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Synopsis: 

The problem of modeling a dynamic system described by a system 
of ordinary differential equations which has unstable components 
Por limited periods of time is discussed. It is shown that the 
global error in a multistep numerical method is the solution to a 
difference equation Initial value problem, and the approximate 
solution is given for several popular multistep integration formulas. 
Inspection of the solution leads to the formulation of four criteria 
for integrators appropriate to unstable problems. A sample problem 
is solved numerically using three popular formulas and two different 
steps! zes to illustrate the appropriateness of the criteria. 



1. Introduction 


When a dynamic physical system is modeled by a system of ordinary 
differential equations 

y' ■ f(y,t) (la) 

y(tQ) - yQ (lb) 

it is possible that the system will be physically unstable for some 

interval of the independent variable t; an example is an aircraft 

approaching a spin configuration. This situation is usually described 

by linearizing the system and saying that it is mathematically unstable 

at t if the Jacobian matrix 3f _ r , ,, 

9y y(t„) y 
1 2 

positive real parts. The usual analysis ’ of a numerical method 
for computing an approximate solution to the system (1) describes the 
behavior of the numerical solution throughout the stable region in 


which fy has all negative eigenvalues. A fundamental result is that 
no multistep formula 


k k 

Ea.y . + hZ 3.f(y .,t .) =» 0 (2) 

. . I n- 1 . - I 'n-i ’ n-i 

1=0 1=0 

for tj = tp + jh, y^ the approximate solution at t ^ , can be stable 
wherever the system (1) is stable and have an error order r greater 
than 2. The error order of a formula is r if the one step truncation 
error is 0(h ), as described in standard texts 

This paper addresses the question of what characteristics a 
multistep formula (2) should have in order to behave well in t-intervals 


2 


when the model equations (1) are unstable. Section 2 gives three 
examples that illustrate the analysis involved in studying typical 
multistep methods. Section 3 states criteria which should be applied 
when choosing a formula to simulate an occasionally unstable system 
and gives an example of an inappropriate choice. Section 4 reports 
a numerical example with three of the formulas analyzed in tiia paper. 
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is in danger of becoming unstable i f hX is large enough and b 0. 

It is easily seen that the leading constant terms 1, -1, are the roots 
2 

of p(z) » z - 1 » 0. When A has a positive real part, is growing 
at the same rate as the correct solution, is being highly damped, 
and the relative error will be almost constant. Thus, the Milne 
method behaves well in unstable regions but is potentially unreliable 
in stable regions of the exact solution. 

A similar analysis can be made of two Adams-Bashforth explicit 
(predictor only) forumulas: 

''n ' Vl t («) 


*“3: y„ - y„., y23f(y„.,,V|) - I6f(y„.2-V2> * 5f(y„.3.t^_3)) (7) 

The error in AB2 is e = + b^^, where e = e , + Xh(3e ,-e -)/2 

n I z n ‘ n"Z 

2 2 1 

which has the characteristic equation ^ - (1 + -^X)C + y hX * 0. By 

assuming distinct roots that are polynomials in hX and equating the 

undetermined coefficients to ( 5 -^ 2 ) > following system of 

equation must be satisfied. If C, ® a. (hX) , and = Z b. (hX) , then 

i-0 ' i=0 ' 

S ^ “ ’’ ^o^ = 0. - 3/2, + b^ = 0, a^b^ + a^b, = 1/2, 

®o^2 ^ Vo ®1^1 = 0, aj + bj - 0, a^b^ + * a,b2 + a2b, = 0, 

which is solved by » 1 + hX + (hX)^/2 - (hX)^^ + O(h^), * 


1 (hX - (hX)^ + (hX)V2) + O(h^). These can be rewritten as 


C, = e 


hX 


(hX)Vl2 + 0(hS, C, - - 1) - 3(hX)2/8 + O(h^). 


5o 


when Re (A) < 0, both error components are decreasing, and ’S 
oscillating. When Re(A) >0, the principal error due to the root' 
of p(C) * 0 is following the true solution, maintaining the same 


relative accuracy. The parasitic error due to the non-unity root of 
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p(C) • 0 is Increasing at a much slower rate due to a zero constant 
term and the exponent hX/2 being smaller than hX. Thus, relative 
error should remain nearly constant in the unstable region of the 
physical system. The relative error will also be good in a restricted 
region of the negative complex half plane called the stability region. 
This consists of all complex hX such that the numerical method produces 
a decreasing solution. If the whole negative half plane is in the 
stability region, the formula is called A-stable.^ Sample stability 
regions for AB2 and AB3 appear in Figure 1. 

The same analysis for AB3 for three separate roots } ( C“C.) and 

i = 1 ' 

for a repeated parasitic root in ) (C"C 2 ) yield contradictions 

in the equations for the undetermined coefficients. Instead of a full 
analysis, then, the worst case can be assumed, that of a multiple 
root corresponding to the multiple root of p(z) = z ( 2 -I) which exists 
for hX =* 0. This means the error looks like + (b+cn) with a 

linearly graving parasitic error term. However, the leading terms of 
^2 be estimated and are found to be ^2 “ ^ ^ O(hX^) = -1) 

+ 0(h ). Thus, even in the worst case of large coefficients a,b,c due 
to large initial errors, the parasitic error is growing at most 
linearly while the solution is growing exponentially, so the largest 
contribution to the error will be the principal error term. 

3. Selection Criteria 

Observation of the three formulas above leads to the following 
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criteria for methods that do well both in the stable and unstable 
regions of the physical model. 

1, All roots of p(z) “ 0 should be very small, and preferably 
zero, except the principal root, which is always 1. In general, a 
nonzero root z. leads to a root of the characteristic equation 

* z.e ; a zero z, yields C- * e - 1, for some real k. 

2, The coefficients k should be positive, so that the parasitic 
error is decreasing for Re(X)<0, and unless the corresponding z. * 0, 

1 k| < 1 is also necessary so the relative error doesn't grow faster 
than the principal error term for Re(X) > 0. 

3, Multiple roots should be avoided, since the error term will 
include factors of n. An exception occurs for z. ® 0, since only 
linear error growth occurs in that case. 

The above three criteria all attempt to keep the error in following 
an exponentially increasing solution from growing faster than the 
solution grows. The following example shows the other extreme. The 
2-step Implicit Backward Differentiation Formula (BDF2) is 

''n ■ J Vi ■ T V2 * f 

and p(z) = 0 has the roots 1, y. The error terms look like 


hX 


+ O(h^), C- * 4 - e + ■|-^(hX)^ + 0(h)^, so this method 


^1 ■ ^2 ’ 3 ■ 27' 

meets criteria 1 and 3, and the small factor 1/3 in lessens the 

problem with the sign of <. However, as can be seen in Fig. 2, the 

stability region is such that the numerical solution, as well as the 

error, is only increasing for a very small portion of the positive 
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half plane. So using BDF2 on an unstable system will likely generate 

a completeiy spurious decreasing numerical solution, except for small 

stepsize h, when the error is increasing fast. A fourth criterion 

k 

has been suggested and will be added here. 

4. The stability region of the formula must not make unstable 
physical solutions numerically stable. 

4. Numerical Examples 

The three methods AB2, AB3, and BDF2 were applied to the unstable 
complex problem 

y' * (1 + j/2)y, (8a) 

y(0) = 1, (8b) 

where j * ("1)^» with two fixed positive stepsizes h and both exact 
and inexact back starting values. The AB2 and BDF2 methods are second 
order methods and thus have comparable errors that are larger than 
those for the third order AB3 method, given the same stepsize. The 
implicit BDF2 method used an Euler predictor and 10 corrector itera- 
tions. The starting values given at n=0 were y^_. * exp(-iXh) for 
the exact values and y^_. = exp(-iXh) + .01 i for inexact starting 
values, where X “ (1 + j/2). All the elements in this example are 
exaggerated so that the results may be easily seen to agree with the 
predictions of the previous sections. The more realistic case of the 
eigenvalues of the Jacobian being close to the imaginary axis, as in 
u spinning aircraft, would require a great deal more computation yet 
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would reach the same conclusions. 

The predicted relative errors can be expressed as 
(a^e + b^\e -1) )e , (a 2 e + (b^ + c^njle -1) )e , 

for AB2, AB3, and B0F2 respectively. 

When h < .25 is chosen, the numerical solution is outside the 

stability region of all three methods. The logarithm of the principal 

.• *i / nhX«nhX\ , • ,, 

relative error is In (a.e /e ) = In a., a constant in all cases. 

If the logarithm of the relative error is not a constant when plotted 

against the independent variable t = nh, then the parasitic errors must 

account for most of this. For the Adams methods and c^ =0, these are 

ln(b. + c.n) (e'^'^^-1 = ln(b. + c.n) + n In(e'^^^-l) -nhX, and 

for BDF2 it is In + n in(y) - j nhX. 

Figure 3 shows the logarithm of the relative error for all three 

methods with exact starting values and h = .05. Note that the BDF2 

error stays less than the AB2 error throughout, indicating only small 

contributions from the parasitic error terms. Figure A shows the same 

results for inexact starting values with a corresponding increase in 

b 2 ,b^ resulting in a large initial error in BDF2. As n increases, the 

term (y)'^e is damped until at t » 4(n*80), this term is again 

hX/2 n 

negligible compared to the b^(e -1) term. Comparison of Figure 3 
and Figure k indicates the c^n (e^ term does not significantly 
affect the error in the AB3 method, except for the initial perturbation, 
although the principal error term is much larger. 

Figure 5 shows the logarithm of the relative error when h * .5, 
which is still outside the stability region of the Adams methods but 



- 8 - 


inside that of BDF2. The error in BDF2 is growing faster than that of 
the other methods since it is due to a decreasing numerical solution 
modelling an increasing analytic solution. Figure 5 comes from exact 
starting conditions but similar rt^sults arise from inexact starting 
values. 

5. Conclusions 

When simulating a problem which is uristable over certain periods 
of the independent variable, one’s choice cf method affects the relative 
accuracy of the numerical solution. When an explicit method is used, 
as is required in real time simulation and is often the case in simple 
simulation programs, a multi step method of the Adams-Bashforth type will 
accurately follow the exact unstable solution despite large initial 
errors. The AB2 method has less computational pitfalls than higher 
order formulas due to the lack of multiple roots, but computational 
experience shows the lack of a leading constant term in these roots causes 
all such methods to be very accurate. Multi step methods are preferable 
to multistage methods, e.g. Runge-Kutta, since a k-step method is as 
accurate with only one function evaluation as a k-stage method requiring 
k evaluations of f(y,t). 

If an implicit corrector is to employed, formulas with a simple 
characteristic polynomial p(z) ■ (l"z)z , i.e. Adams-Moul ton 
methods, should behave better than even the (possibly) A-stable BDF 
methods because they have only one non-zero root and produce increasing 
numerical sequences wherever the solution is increasing, at least for 
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methods of less than five steps. The best implicit method from a 
stabi I i ty/instabil i ty standpoint Is the trapezoidal formula 

''n ■ Vl * 

which has no parasitic roots and is increasing (decreasing) only 
where the exact solution is unstable (stable). However, this method 
is only second order. All of the implicit methods take on their opti- 
mal stability properties only when iterated to convergence, which 
usually requires at least two evaluat ions of f(y,t). Therefore, from 
the standpoint of efficiency and stability, AB2 seems to be a practical 
integration method for occasionally unstable systems of ODE's. 


6. Acknowledgements. 

This work was supported in part by the National Aeronautics and 
Space Administration through its grant NSG 1335. 


References 

1. Henrici, P., 'Discrete Variable Methods for Ordinary Differential 
Equations,* Wiley, New York, 1963. 

2. Gear, C. W. 'Numerical Initial Value Problems in Ordinary Differ- 
ential Equations,' Prentice-Hall, Englewood Cliffs, N. J., 1971. 

3. Dahlquist, G., "A Special Stability Problem for Linear Multistep 
Methods," BIT, 1963 , 3, PP. 27-^3. 

4. Lindberg, B., "On a dangerous Property for Methods for Stiff 
Differential Equations," BIT, 197^, pp. ^30-436. 


Captions to figures 


1. Stability regions for AB2 (larger figure) and AB3. 

2. Stability region for B0F2 (method is unstable only inside the 
region indicated). 

3. Logarithm of relative error for h « .05, exact starting values - 
AB2 (12), AB3 (A), B0F2 ({>). 

4. Logarithm of relative error for n ^ .05, inexact starting values - 
AB2 (□), AB3 (A), BDF2 (^). 

5. Logarithm of relative error for h ■ .5 " AB2 (□), AB3 (A), 

BDF2 {(}), 
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